o 

(N 



00 

in 

o 



X 



Noname manuscript No. 

(will be inserted by the editor) 



Sequential Convex Programming Methods for Solving 
Nonlinear Optimization Problems with DC constraints 



TranDinhQuoc Moritz Diehl 



Received; date / Accepted: date 



Abstract This paper investigates the relation between sequential convex programming (SCP) 
as, e.g., defined in |24| and DC (difference of two convex functions) programming. We first 
present an SCP algorithm for solving nonlinear optimization problems with DC constraints 
QQ ■ and prove its convergence. Then we combine the proposed algorithm with a relaxation tech- 

P\| I nique to handle inconsistent linearizations. Numerical tests are performed to investigate the 

behaviour of the class of algorithms. 

\^ ' Keywords Sequential convex programming ■ DC constraint • relaxation technique • 

(7^ , nonconvex optimization. 

1 Introduction 

Let io(R") denote the set of all proper lower semi-continuous convex functions from R" to 
R, and ^'^(R") := ro(R") -ro(R") denote the set of DC functions on R". We are interested 
in the following nonconvex optimization problem: 

min fix] 
xeR" ^ ' 



s.t. g{x) < 0, (P) 

xeQ, 

where / : R" — > R is convex, Q is a nonempty closed convex subset in R", and ^ : R" — > R"' 
withg = (gi,...,gm)^ andgi (;' = 1, ... ,m) belongs to ^"^(R"). We refer to g(jc) < as DC 
constraints. Let us denote by D := {x € Q : g{x) < 0} the feasible set of ^ and intD the 
set of interior points of D. 

Problems of the form ^ have been studied by many researchers in theory and applica- 
JH , tions (see, e.g., I2 HTT1I 1211 1311 1411291 and the references quoted therein). However, the meth- 

ods for solving (O that exploit DC structures are usually global optimization techniques. 
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These approaches are not applicable to problems with a high dimension. In this paper, we 
are interested in finding local minimizers only. 

The class of DC functions is sufficiently rich to deal with many practical problems. It is 
well-known II13II14I that the set of DC functions defined on a compact convex set of R" is 
dense in the set of continuous functions on this set. Therefore, in principle, every continuous 
function can be approximated by a DC function with any desired precision. Moreover, every 
C^-function defined on a compact set is a DC function 1 1 1 1 that includes the smooth cases 
of ^. Many practical problems can be reformulated in the form of ^ (see, e.g., 1 1_4J). 
Although DC representations are available for important function classes, finding such a 
representation for an arbitrary DC function is still a hard problem. 

This paper investigates the relation between SCP methods I18II24II and DC programming 
IT]|2l l23l . Both families of methods address the local solution of nonconvex optimization 
problems via an iteration based on convex subproblems. 



1 . 1 DC programming 

DC programming algorithms (DCA) for solving l|Pll have been introduced by Pham (Tl|2l 
|23l. The original DCA is supposed to solve convex constrained DC programs. To handle DC 
constraints, penalty functions have been used L21 and then DCA is applied to the penalized 
problem for a fixed penalty parameter. Yuillie and Rangarajan in 1301 proposed a method for 
solving smooth DC programs that is called the concave-convex procedure (CCCP), a variant 
of DCA applied to smooth DC programs |27|. The authors in |28| further investigated the 
global convergence of the CCCP method. DCA as well as CCCP have been widely applied 
in many practical problems (see, e.g. ITl l27ll30l ). It is well-known that the use of penalty 
functions in DC programming with DC constraints introduces conservatism and might lead 
to excessively short steps. 

One particular variant of DC programming that keeps the DC constraints in the problem 
was considered in 1261 . This again leads to possible conservatism or even to infeasibility of 
the subproblems (which might be overcome by relaxation techniques). These methods have 
not become very popular due to these problems and their combination with exact penalties 
was never properly investigated. 

It is the aim of this paper to improve and investigate the numerical behaviour of these 
algorithms and show that they can be interpreted as a special case of SCP methods. 



1.2 Sequential Convex Programming 

In 1241 . a generic algorithm framework for solving nonlinear optimization problems with 
partially convex structure was proposed which is called sequential convex programming 
(SCP). The main idea of SCP methods is to convexify the nonconvex part and preserve the 
remaining convexity in the resulting subproblems at each iteration. Under mild assumptions, 
the local convergence of the SCP methods was proved. The rate of local convergence is 
linear. 

To the family of SCP methods belong such classical algorithms as the constrained or 
unconstrained Gauss-Newton methods as well as sequential linear programming (SLP) or 
sequential quadratic programming (SQP) with convex subproblems I10lll5ll20l . All these 
methods are based on linearization of nonconvex constraints or objective functions, and are 



widely used in applications of nonlinear optimization, in particular, in parameter estimation 
(constrained Gauss-Newton |3| and nonlinear model predictive control |8,£]). 

When DC constraints are treated within an SCP framework, it is possible to only lin- 
earize the concave parts. This can be interpreted as a special case of SCP, which offers 
a favourable feature: namely that globalization strategies like line search or trust-region 
methods are not needed and full SCP steps can always be taken. When feasibility of the 
subproblems becomes an issue, which is always the case for nonlinear equality constraints, 
we propose to relax the subproblems using an exact Li -penalty function and investigate the 
behaviour of this relaxed SCP algorithm. We show through an example that it can lead to 
less conservative convex subproblems than the standard approach of using unconstrained 
DC programming with penalty functions. 



1.3 Notation and definitions 

Throughout this paper, we use R'" for the set of m-nonnegative vectors and R+ (resp., R^) 
for the set of nonnegative (resp., positive) numbers. 

A function / : R" — )• R is called p-^ -convex on a convex subset X of R" with p^ 6 R+ if 

for all x,y eXandte [0, 1] the inequahty f{tx+{\ - t)y) < tf{x) + {l-t)f{y) - ^f (1 - 
?)||jc — y||^ holds. If p-'^ = then / is convex. Otherwise, / is strongly convex with the 
parameter p^ > 0. 

Let us assume that / is a DC function such that f = f\ — fi, then it is trivial to see 
that / = (/i -K £|| • f) _ (/2 + £|| . ||2) foj- any given p > 0. Therefore, without loss of 
generality, we can find a DC decomposition (/i,/2) of / such that f\ and /2 are strongly 
convex. We also use the notation dom/ := {x 6 X | f{x) < +o°} for the domain of a convex 
function /. For x E dom/, the symbol df(x) denotes the exact subdifferential of / at x, i.e., 
df{x) := {£, 6 R" I f{y) > f(x) -(- (^^(y - x), \/y E X}. A convex function / is said to be 
subdifferentiable at jc e dom/ if df(x) / 0. A vector E, 6 df{x) is called a subgradient of / 
atx. 



1 .4 Optimality condition 

Suppose that (m, v) is an arbitrary DC decomposition of g. Let us define L{x, A ) := Ao/(x) -|- 
A^[m(jc) — v{x)] the Lagrange function of problem (|P]|. The generalized F. John condition of 
^ is expressed as follows (5): 



£ Ao<9/W + I^-li ?ii[dui{x) - dvi{x)] +Nn (x), 
0^(Ao,A)>0, u{x)-v{x)<0, A^[m(x)-v(a:)]=0, 



T, „ (1) 



where df{x), dui(x) and d\'i(x) (i = 1, . . . ,m) are the subdifferentials of /, m,- and v,- at x, 
respectively. The multivalued mapping A^^^ is the normal cone of 12 at x defined by: 

[0 otherwise. 

Note that the first line of lUll includes implicitly that x 6 12 . If (x* , Aq , A * ) satisfies (O then 
X* is called a stationary point and (Aq , A*) is the corresponding multiplier of ©. 



Since problem © is nonconvex, a stationary point might not be a local minimizer. How- 
ever, we will show later that under the calmness constraint qualification, the first order nec- 
essary condition for l|Pll still holds. 

We consider the following parametric optimization problem: 

V{5) := M{f{x) I g{x) <5,xen}, (P(5)) 

where the perturbation (or parameter) 5 belongs to a neighborhood Ue C R™ of the origin. It 
is trivial that P(0) = P. Let x* solve (O. Problem ^ is said to be calm at x* (in the sense of 
Clarke's calmness constraint qualification [5]) if there exist a neighborhood Ue of the origin, 
Xe ofx* and a positive number T such that for all 5 e [/^ and xeXe that are feasible for |P(3)[ 
one has f{x) — f(x*) + t||5|I > 0. The characterizations of calmness have been investigated 
in the literature (see, e.g., 1 5 , 16,^|). The optimality conditions for DC programs with DC 
constraints have been studied in ITtI . 

If v/ (;' = 1, . . . ,m) is continuously differentiable on R" then, under the calrrmess of (|P]| 
at a local solution x*, without loss of generality, we can assume that the multiplier Ao = 1. 
Thus the F. John condition lO collapses to the (generalized) KKT condition. With Aq = 1, the 
point (x*, A*) satisfying ([TJ is called a KKT point. In particular, if/, m; and v,- (/ = 1, . . . ,m) 
are continuously differentiable on R", and Q is the whole space, then the condition ([TJ 
collapses to the classical KKT condition in smooth nonlinear optimization |20]. Under the 
Mangasarian-Fromowitz constraint qualification, the first order necessary condition corre- 
sponding to (TJ holds for (|Pj. The following theorem shows that the first order necessary 
condition for problem (|Pj still holds. 

Theorem 1 Suppose that f £ F(R") and («,v) is a DC decomposition of g such that v is 
continuously differentiable on R". Let x* be a local minimum of (|Pj such that (|Pj is calm at 
x*. Then there exists a multiplier A* £ R"" such that {x* ,X*) is a solution to the KKT system 

Q. 

Proof Note that if a function (p is continuously differentiable (resp., convex) then the Clarke 
subdifferential coincides with its gradient (resp., its convex subdifferential) 1 5 1 [Proposition 
2.2.7]. Since v,(-) is convex and continuously differentiable on R", it implies that dvi = 
{Vv,} for all i = 1, . . . ,m. On the other hand, since m, is subdifferentiable on R", we have 
d'^gi = <5'^("( — V,) = dui + V(— V,) = dui — Vvj, where d'^gi is the Clarke subdifferential 
of g; (i = 1, . . . ,m) Is]]- Applying Proposition 6.4.4 in (5) we obtain the conclusion of the 
theorem. D 

The rest of the paper is organized as follows. Section |2] presents two motivating exam- 
ples. A variant of the SCP algorithm for solving (|Pj is presented in Section[3] Then its global 
convergence is investigated in Section|4l A relaxation technique is proposed in Section|5]to 
handle possibly inconsistent linearizations. Computational tests are performed in the last 
section to demonstrate the behaviour of the class of algorithms. 



2 Motivating examples 

There are many practical problems that can be conveniently reformulated in the form of (|Pj 
such as mathematical programs with complementarity constraints, bridge location problems, 
design centering problems, location problems, packing problems, optimization over efficient 
sets, trust-region subproblems in SQP algorithms, and nonconvex quadratically constrained 



quadratic programming problems (see, e.g., I13II23I ). For motivation, we present here two 
examples. The first example originates from optimal control of a bilinear system and the 
second one is a mathematical programming problem with complementarity constraints. 



2. 1 Nonlinear model predictive control (NMPC) of a bilinear system 

The optimization problem resulting from NMPC of a bilinear dynamic system has the fol- 
lowing form: 

niin Fo(x, u) := i if^ ' [x[W*x, + M[lV>i] + ^xJ^W^xh^ 
s.t. x/c+i =Axk+B[xk,uii] +Cuii, k = 0,...,Hp — l, 

-'^0=JCinit, _ (NMPC) 

Xj,<xi,<Xk, k = 0,...,Hp, 
M.k ^ '^k ^ Uk, k = 0,. . . ,Hp — I, 

x]jWeXHp < rf, 

where W^, W*, We are the weighting matrices; A,C are given consistent matrices; Xi„it is 
a given initial state; Xj^, xi^, Uj^, Uj^ are lower and upper bounds on the variables x^. and m^., 
respectively; r/ > is the radius of the terminal region; and B[xk , Uk] denotes a bilinear form 
of ;<:«: and u^. 

Introducing a new variable w : = (xg ,xj ,. ..,xjj , Mq , . . . , mJ - 1 ) ^ ^ ^"" ^i^h n^f = {Hj, + 
Vjiix + HpHii, the objective function of ( INMPCb can be rewritten as Fq{w) = jw^Z/w, where 
// is a symmetric positive semidefinite matrix with W^ , W* and W^ on the diagonal block. 
It is known that a given bilinear form is always associated with a quadratic form. Therefore, 
the discrete time bilinear dynamic system x^^+i = Axk +B[xk,Uk] +Cuk (k = 0,. ..,Hp— I) 
can be reformulated as: 

w^PiW + qfw + ri = 0, i=l,...,m, (3) 

where m := HpH^, P, is a given symmetric indefinite matrix, qi e R"" and r,- e R (; = 
l,...,m). Any symmetric indefinite matrix Pj can be decomposed in such a form P, := 
P^ — Pf, where Pf and Pf are two symmetric positive semidefinite matrices (e.g., using 
spectral decomposition). Using two different DC decompositions of Pi and choosing q], qj, 
q\ , §?, r] , rf, f\ , ff such that qi = q] — qf = q\ — qf, r,- = rj — rf = r] — rf, respectively, the 
equality constraints ([3} can be rewritten as 

[{w^P^w+ {q]fw + rj]- [w'^Pfw + {qffw + rf] < 0, 
[{w'^Plw+{qffw + ff]-[w'^Pfw+{qjfw + rj] < 0, 

for all i = 1, . . . ,m. Hence, problem I INMPCb is reformulated in the form of ^. Note that 
Pj = pf (j = 1, 2) is a possible choice in the formula ^. 



2.2 Mathematical programs with complementarity constraints 

Mathematical programs with equilibrium constraints (MPEC) have been studied widely and 
have many applications in economic models, shape optimization, transportation, network 



design, and data mining. In this example, we particularly consider the following mathemat- 
ical programming problem with complementary constraints: 

min f{x,y) 

s.t. {x,y) e S, (MPCC) 

x>0, Cx + Dy + e>0, 
x^{Cx + Dy + e)=0, 

where / : R" x R'" — > R is convex, S C R"+'« is a nonempty closed convex set, e 6 R", and 
C, D are two given matrices of consistent dimensions. 

Theory and methods for dMPCCI l have been developed intensively in recent years (see, 
e.g., ||6| |19||221 and the references quoted therein). The main difficulty of this problem is 
the complementarity constraints in the two last lines of JMPCCt . These constraints lead to 
nonconvexity and loss of constraint qualification of the problem. 

Introducing a slack variable z, the complementarity constraint can be reformulated as: 

Cx + Dy-z + e = 0,x>0, z>0, x'^z = 0. (5) 

Since x>0 and z > 0, the constraint x^z = is equivalent to x^z < 0. Using the expression 
2x^z = \\x\\^ + \\z\\^ — \\x — z\\^, we can rewrite the condition x'^z < as a DC constraint: 

u{x,z)—v{x,z)<0, (6) 

where u{x,z) := |l(^,z)|l^ and v(x,z) '■= ||j: — z|p that are convex. Problem ( IMPCCb is now 
reformulated in the form of ^. 

For an MPEC problem, by using the KKT condition for the equilibrium constraint (low 
level problem), we can transform this problem to the form (MPCC) (see [6]). Then, by the 
same technique as before, we obtain a DC formulation for the equilibrium constraint. 



3 Sequential convex programming algorithm with DC constraints 

In this section, we present an algorithm for solving problem ^ which we might call se- 
quential convex programming with DC constraints. Let us assume that (m,v) is a DC de- 
composition of g, i.e., 

g{x) = u{x)-v{x). (7) 

For a given point x* 6 12, we take an arbitrary matrix S* 6 dv{}f), where the multivalued 
mapping dv{J^) := [dv\[J^Y , . . . ,dv,„[x''YY with dvi[J^) (i = l,...,m) is the subdiffer- 
ential of the convex function v,- at x*. We will refer to S*^ as a subgradient matrix of v at jc*. 
Consider the following convex problem: 

min fix) 

A-eR" ■' ^ ' , 

s.t. m(x)-v(x*)-S*(jc-/)<0, (P(-«)) 

xe Q. 

Since |P(x^ is convex, under the Slater constraint qualification 

n{n)nlx:u{x)-v{J')-E''{x-x'')<o\ ^0, (8) 



w here ri jQ ) is the set of relative interior p oints o f the convex set Q , any global solution x^^^ 
of |P(x^)| is characterize d as a K KT point of P(.t^') In the following algorithm, we assume that 



the convex subproblem P(;c*) is solvable for given x'' and S*. 

A generic framework of the sequential convex programming algorithm with DC con- 
straints (SCP-DC) can be described as follows: 

Algorithm 1 

Initialization: Take an initial point jc*^ in 12. Set k := 0. 
Iteration k: For a given J', execute the three steps below: 

Step 1: Compute a subgradient matri x 5^ £ (9v(x*) 



Step 2: Solve the convex subproblem P(x*) to get a solution .«'^+' and the corresponding 
multiplier A'^^'. 

Step 3: If ||;c'^+' — Jf*|| < e with a given tolerance e > 0, then stop. Otherwise, increase k 
by 1 and go back to Step 1. 

At Step 1 of Algorithm[T] a subgradient matrix S*^ of v at x* must be computed. If v,- (/ = 
1, . . . ,m) has a simple form, S* can be computed explicitly. Otherwise, a convex problem 
needs to be solved. If v is differentiable at jc* then dv{jJ') is identical to the Jacobian matrix 
of V at jc*, i.e. dv{x'') = {Vv(x*)}. 

The cost of finding an initial point jc" 6 12 depends on the structure of 12 . It can be 
computed explicitly if Q is simple. Otherwise, a convex problem should be solved. The 
projection methods (onto Q) can be also used in this case. 

Remark 1 If the objective function / of (IQ is linear (resp., quadratic) then: 



i) If the function u is linear then subpr oblem P(x ) is linear (resp., quadratic) 



ii) If the function u is quadratic then |P(jc*)J is a quadratically constrained linear (resp., 
quadratic) programming problem. This problem can be reformulated as a second order 
cone programming or semidefinite programming problem f4\. 

DC decomposition of the function g plays a crucia l role in Algorithm [T] A suitable DC 
deco mpositi on may ensure that the convex subproblem |P(;c*)| is solvable. Moreover, it might 
make P(x*) easy to solve, and help Algorithm [T] to reach a KKT point of (O (e.g., u and v 
have small strongly convex parameters). The following small example shows the behaviour 
of Algorithm[T]using two different DC decompositions. 





min f{x) := —Ax\ -\-X2 
' s.t. g{x) ■.= x\-xl-A<Q, 


(9) 




y xeI2:=[-3,3]x[-2,2]. 




The constraint x\—x\- 


- 4 < is a DC constraint. Hence, for 


a given tolerance e = 10"^ 


and a starting point x^ - 


= (0,0)^, if we choose (m, v) with u{x) 


= x\ and v{x) := x\ for the 



DC decomposition of g (Case I) then Algorithm [T] converges to the global solution after 2 
iterations. If we choose u{x) := x\ +x\ and v{x) := 2x\ (Case 2) then it converges to the 
global solution after 4 iterations. Note that, in the first case, u and v are only convex, while 
u is strongly convex with the parameter p" = 2 and v is convex in the second case. The 
convergence behaviour is illustrated in Figure [T] Here, the left figure corresponds to Case 1 
and the right one corresponds to Case 2. 

The following lemma shows that if Algorithm[T]terminates after some iterations then x* 
is a stationary point of (JFj. 






Fig. 1 Convergence behaviour of Algorithm^using different DC decompositions. 



Lemma 1 Suppose that x is a solution o/IPfx^jf then it is a stationary point of the original 
problem dPl l. 



Proof Suppose that }f is a solution of P(j(^) corresponding to the multipHer A*^ then (;i;*, A*^) 
is a solution of its KKT system, i.e., x* e\o, 6 df{)f) + [du{x^) - S*^]^A* + A?i2(jc*), 
m(/) - v(/) - (S*)(/ -/) < 0, A* > and (m(/) - v(jc*) - E'')(x^ -x^)fX'' = 0, which 
implies that jc* e 12, e df{x^) + [<9m(/) - (9v(x*)]^A* +A^i2 (^*), m(^*) - v(;t:'^) < 0, A* > 
and [m(jc*) — v(jc*)]-^A'^ = 0. The last five relations mean that {x'',?i'') satisfies (|T]|. Thus a^ is 
a stationary point of ^ corresponding to the multiplier A*. D 



4 Global convergence of the SCP algorithm with DC constraints 

The next lemma gives us a key property to prove the global convergence of Algorithm [T] 

Lemma 2 Suppose that f, m; and v; (/ = 1 , . . . , m) are p^ , p"' and p '' - convex, respectively. 
Then the sequence {(;c , A )} generated by Algorithm\l\satisfies 



/(/)-/(/+') >l(p/ + J^p".A/+')||/+'-/| 



+^Ep"'AMx^-/ 



(10) 



-1||2 



2' 



/=i 



Proof S ince jc^ +^ is a solution of P(x'^) corresponding to the multiplier A*+', the KKT con- 
dition of ] 



P(x ) is expressed as follows: 



> m(jc*+1) - v(x*) -S*(;c*^+i -x^), A'^+l > 0, 
= (A*+')^[m(/+1) - v(/) - E\x^+^ -/)]. 



(11) 



From the first line of dllb . we have 

(^)+')^(y-/+') + (A*+')^[S,fi -S^](y-/+i)>0,VyeI2, (12) 



for all vectors (§J:+' 6 (9/(x*+') and matrices S*+' e du{)^ 



+h 



Since / and w/ (/ = 1 , ■ • • , fn) are strongly convex on ^ , it holds that 

/(3')-/(^'+')>(^f')^(y-^'+') + Ylb-^+'ll',Vyei2, (13) 

M(y) -«(/+') >S,^'(y-/+i) + ^|b-/+'|j2,Vy 6 12, (14) 

where p" = (p"i , . . .,p""'Y . Combining ( fT2] |. ( fT3] | and ( fT4l l. and noting that A*+' > 0, we 
obtain 



1 m 

/—I 



^ , w „ ,mV^+^ — -^11^ 



2^ ,=1 



(15) 



1 Hi 

^ /-I 

1 «7 

>2[P^ + EP"'A'+']lb--^+'f,V}'eI2. 

Substituting _y = x^ 6 ^ into ( fTSl l and after a simple rearrangement, we get 

/(/) + (A^+')^[«(/)-v(/)] 

-/(-«^+') - (A'=+')^[m(/+') - v(/) - S'=(/+i -/)] (16) 



Using the third line of jllb . the inequality ( 116b is reduced to 

1 m 

/(/) + (A^+')^[m(/)-v(/)] -/(/+') >^[p/ + £p"'Af+i]||/+i-/f. (17) 



Now, since v/ (;' = 1, . . . ,m) is p'"' - convex, we have 

v{/+')-v{j^)>E\/+'-x') + ^\\^+'-x'\\\ 
where p'' = (p ''',.. . ,p^''"Y . This inequality implies that 

m(/+1) - v(/+i) < m(/+1) - v(/) - S*(/+i -/) - ^||/+' -/||^ (18) 

Using the second line of jllb for jlSI l, we obtain 

M(/+')-v(/+i)<-^||/+'-/f<0. (19) 

Applying ( |19l l with x*^ instead of x*+' to ( 117b yields 

1 m 1 m 

/(/)-/(/+') >^[p/ + Ep"'A'+']ll-^+'--^ll' + ^Ep"'A'+'ll^-^'"'ll', (20) 

which proves l llOb . D 
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Remark 2 From the proof of Lemma [2] (see l ll9b ) we can see that Algorithm [T] always gen- 
erates a feasible sequence {x''} to (O. If p^' > then it is strictly feasible. Thus AIgorithm[T] 
can be considered as an inner approximation method. 

Remark 3 If either / is strongly convex or at least one function m,- (reps., v,) {i = 1, . . . ,m) 
with respect to A,. ' > is strongly convex then the sequence of the objective values {/(jc* ) } 
is decreasing. 

The convergence of Algorithm[T]is stated by the following result. 

Theorem 2 Suppose that f is bounded from below on D, and the sequence {(jc , A )} gen- 
erated by Algorithm\^is bounded on £2 x R™. Then: 

(i) Ifp' > thenlimk^oa ^x — x^|| =0, and every accumulation point (x*,A*) o/{(jr, A )} 

is a KKT point of (O. 
(ii) If there exists an index (q G {1, • • • ,m} such that p"'o > (resp., p''(i > 0) then 

limAf+i||/+i-/||2=0 {resp., lim 1/+' ||/ -/-' f = 0), 
and every accumulation point o/(x*,A*) of {{x ,?i )} such that X^ > is a KKT point 

ofB- 

(Hi) If the set of the KKT points of ^ is finite then the whole sequence {(jr,A )} converges 
to a KKT point of ©. 

Proof From Lemma|2l it turns out that the sequence {f{x^)} is nonincreasing and is bounded 
from below by assumption. Then it converges to /* > — o°. Summing up inequality dlOb from 
k= I to k = N and then passing to the limit as fc — > o° we obtain 



I 



1 m 1 «; 

^ !=1 ^ j=l 



<f{A-r<+° 



(21) 

If p-^ > then the inequahty dlTT l implies that lim^.^^ ||jc'^+' — x*|| = 0. Since {(jc*,A*)} is 
bounded by assumption, it has at least one limit point. Suppose that (x*, A*) is a limit point 
of {(x*,A*)}, which means that there exists a subsequence {(jc*,A'^)}^.gj^ of {(jr*, A*)} such 
that (x*, A*)(/t e JT) ^ (x*,A*), where A* 6 R"^. Since df, dui and dv, (i = 1, . . . ,m) are 
upper semicontinuous, passing to the limit of the subsequence as A:(e J(f) — )• oo in dllb we 
conclude that (x*, A*) is a KKT point of (jQ. The statement (i) is proven. 

For the statement (ii), it is sufficient to prove the first case (i.e., there exists /q such that 
p"'o > 0), the second case is done similarly. Suppose that there exists at least one index (q G 
{1, . . . ,m} such that p"'o > 0. Using again ( l2Tb . it is easy to show that lim^^„ A/"'"' ||jc'^+' — 
x*|p = 0. As before, if (x*,A*) is a limit point of a subsequence {{J^,^'^)}ke,je such that 
A* > then we have lim^.(g^Wt„ 11^^' ~-*'^ll = 0- Passing to the limit through the subse- 
quence as fe(6 J(f) — 5- o° in illlt we conclude again that (x*, A*) is a KKT point of ^. 

The last statement (iii) can be proved similarly using the same technique as in l21i rChapt. 
28]. D 

Suppose that x* is a stationary point of (O associated with a multiplier A*. If we denote 

by 

/+(x*):={(e{l,...,m}|A*>0} (22) 

the strictly active set of (O at x* , then the assumption (ii) in Theorem[2]requires that 7+ ix* ) ^ 
and at least one function m,- (or v,) / £ I+{x*) is strongly convex. 
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Remark 4 {Regularization). From Lemma[2l we see that if/, m,- and v,- are only convex for all 
i £ I+(x*) (but not strongly convex) then Algorithm[T]might not make / strictly decreasing, 
i.e, f(pf^^) / f{^) for fc > 0. In order to overcome this issue, a regu larizat ion term can be 
added to the objective function of P(x*) Instead of solving problem |P(-y*)| Algorithm [T] is 



modified at Step 2 by solving the following regularized problem: 

f min/(x) + £||;c-/|P 

S.t. u{x)-v{J')-E\x-x'')<0, (PrK)) 

where S* 6 (9v(x*) is arbitrary, and p > is a regularization parameter. This technique is 
closely related to the proximal point methods 11811251 . 

However, using the regularization term with a large p may lead to short steps. Consequently, 
Algorithm [T] converg es slow ly to a KKT point. In practice, we only add the regularization 
term if the solution of 'P{]f) does not make / strictly decreasing at the current iteration. Note 
that if p > then Algorithm [T] always makes / strictly decreasing, i.e., /(x*) — f{J'^^) < 
-f||/+i-.=c*f<Ofor;c^+i^jc*. 

Remark 5 {Handling the DC objective function). If the objective function / of (IQ is also a 
DC function and f[x) = /i (;i:) — /2(x) is a DC decomposition of /, then subproblem |P(j:* )| 
at Step 2 of Algorithm[T]is replaced by the following convex subproblem: 



min 

V6R" 


/iW- 


-fii^) 


-(^A 


Y{x 


-/ 


s.t. 


u{x)- 

xen 


v(^)- 


S*(x 


-x^) 


<0 



(Pdc(x*)) 



with matrix S* e dv{J') and vector <§i £ (9/2(x*). The conclusions of Theorem |2] are still 
valid for this modification. A smooth variant of this algorithm was considered in 1261 applied 
to DC programs arising in support vector machines, without convergence theory, however. 



5 A relaxed SCP algorithm with DC constraints 

According to DCA approaches, to handle DC constraints, a penalty function is used to bring 
these constraints into the objective function ]2A- The obtained problem becomes an uncon- 
strained or convex constrained DC program, and the unconstrained DCA can be applied to 
solve this problem. We start this section by introducing one possible DC decomposition to 
handle the DC constraints using Lj -penalty functions, which is often used in practice | OH- 
We will show through an example that by using an Li -penalty function to handle the DC 
constraints, DCA may make only slow progress to a stationary point of ©. 
Let us define the Li -penalty function of ^ as follows: 

^{x■^l)■.= f{x)+^i\Mx)]+\\^, (23) 

where /i > is a penalty parameter and [g(ji:)]+ = max{g(x),0}. Note that if g has a DC 
decomposition (m, v) then we have [g(-«)]+ = vt\ax{u{x) ,v{x)} — v{x). Since u and v are con- 
vex, max{M,v} is also convex. Thus (max{M(-),v(-)}, v(-)) is a DC decomposition of [g(-)]+. 
Since(^(x;/x)=/(j:)-F/xI"li[max{M;(x),v,-(x)}-v;(x)]=/(x)-F/iI™,[max{M,-(x),v;(.ic)}- 
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/X Vj (x) ] , if we define Ufj{x) : = / (x) + /i Yi'lL i max { m; (jc) , v; (j:) } and v^ (x) : = jU YUL i v; (x) then 
(^ (x; jU ) is a DC function, and (m^ , v^ ) is a DC decomposition of (j){x;^). 

Tlie Li -penalized problem associated with (O can be rewritten as a convex constrained 
DC program: 

mm{(l>{x;n)=u^{x)-v^{x)}. (?);'=) 

DC A f2l starts from an initial point x" e 12 and generates a sequence {x*} by solving the 
following convex subproblem: 



mmM 

xen 



„(x)-v^(/)-(^*)^(x-/), (P-(x*)) 



where S,^ 6 (9v^ (x*) and /i is fixed to a suitable large value. It is proved in (2) that for this 
DC decomposition, there exists an exact penalty parameter /i/ > such that for all /i > /ii. 



any solution of problem 1?^"^ I solves ^. 

Now, we show that by using this standard technique, DCA may lead to slow convergence 
to a stationary point. Indeed, we consider an example by minimizing a convex function / 
subject to a DC quadratic constraint ^{x^Px — x^Qx) < p^x + r, where matrix P is sym- 
metric positive semidefinite, Q is symmetric positive definite, p £ R", and r £ R. If we 
define m(x) := jx^Px — p^ x — r and v(x) := ^x^ Qx then v is strongly convex with parame- 
ter p'' = A]nin(2), where XminiQ) is the smallest eigenvalue of Q. Applying DCA to problem 



(P"'^ I we have v^(x) = piv{x) that is strongly convex with parameter p''^ = /jAmin(G)- If M 
IS large then p'*' is also large. In this case, DCA makes only slow progress to a stationary 
point of ((Pji. 

Instead of using the penalty function i l23t directly, in the SCP framework, we automati- 



cally obtain a different relaxed algorithm. We first relax the DC constraints by 

m 

\mB.f{x)+^Ysi 
(—1 

s.t. m(x)-v(x) <s, (24) 

xeI2,i>0. (25) 

We use a relaxation technique to handl e possi bly inconsistent linearizations that may lead to 



infeasibility of the convex subproblem P(x*) in Algorithm|T] Note that m(x) — s is convex in 



{x,s) as well as v(x). Each SCP-DC subproblem is then given by: 

min/(x)+Mirii^/ 

s.t. m(x)-v(x*)-S*(x-x*)<.?, (P(j^;M)) 

.5 > 0, X 6 12 . 

A relaxed variant of Algorithm[T]called relaxed SCP algorithm with DC constraints (rSCP- 
DC) is described as follows: 

Algorithm 2 

Initialization: Choose a penalty parameter fio > 0. Take an initial point x" in 12 . Set k :=0. 
Iteration k: For a given x*, execute the three steps below: 

Step I: Compute a subgradient matri x 5^ £ dv {x 



Step 2: Solve the convex subproblem P(x*; /i) with jU = /i^ to get a solution (x*^+' , i*+^ ) 
and the corresponding multiplier A*+T 
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Step 3: If ||j(^+' — x!^\\ < e and ||i*+'|j < e with a given tolerance e > 0, then stop. 
Otherwise, update the parameter fi^, increase A: by 1 and go back to Step 1. 



Note that the subproblem P{x'';^) is always feasible and the convergence theory of the pre- 
vious section is applicable. However, the parameter /x influences the behaviour of Algorithm 
|2] If the parameter fi is chosen too large, the minimizati on enforces s to decrease, which re- 



duces the infeasibility gap of the subproblems P(x ,jU) too fast. Otherwise, the infeasibility 



gap s may be increased. Balancing between the optimality and the infeasibility plays an 
important role in Algorithmic The parameter fit can be fixed to a "suitable" value or up- 
dated at each iteration of the algorithm. A refined variant, which is however not the topic of 
this paper, separately updates penalty parameters /i,- for each j, and make sure that they are 
sufficiently large, but not much larger than the corresponding constraint multipliers. 

The following inequality shows that Algorithm |2] makes a decreasing progress of the 
objective function /^ {x, s) := f{x)+fj. Y!ILi ■s/- 

Corollary 1 Suppose that f, ui and v,- (;' = !,...,/«) are p' , p" and pj - convex, respec- 
tively. Then the sequence {(x ,A ,.s )} generated by Algorithm\2\satisfies 



1 "' 

/^(/,/)-/^(x^+i,/+') > (p/ + £p«.A/+')||/+i-/|, 

^ (=1 



2 



1 m 



(26) 



2,=i 



where f^{x,s):=f{x)+fil,1L^ Si. 



The conclusions of Theorem [2] still hold for this case, where the objective function is 
ffi (x, s) (with a fixed value fi > 0) instead of /. 



6 Numerical tests 

To verify the performance of Algorithms [T] and |2] we implement two numerical examples. 
The first example solves nonconvex quadratically constrained quadratic programs (ncvQCQP). 
The second one is a mathematical program with complementarity constraints. 



6.1 Example 1 

Consider the following indefinite quadratically constrained quadratic programming prob- 
lem: 

min fix) := Ix^Qx + q^x 

xeR" •' "■ ' -^ 

s.t. \x^Px + p'^x<a, (ncvQCQP) 

Ax <b, 
[<x<u, 

where q,p,l,u eW, a 6 R, ^ 6 R'"^, g is a symmetric positive semidefinite matrix in 
R"^", A is an m2 x n real matrix, and P is an « x « symmetric indefinite matrix. If P is sym- 
metric positive semidefinite then problem l |ncvQCQ"Pl l is a convex quadratically constrained 
quadratic programming problem (QCQP) I4j. 
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Table 1 Computational results of Algorithms ^andlJlfor jncvQCQPj. 



Problem Information 


Algorithm LJ 


Algorithmic) 


FMINCON 1 


ni2 


n 


/* 


error 


iter 


time 


iter 


time 


M* 


p 


iter 


time 


5 


10 


121.3768 


3x10-^* 


64 


16.13 


3 


1.88 


0.1 





36 


1.54 


10 


30 


12.1228 


3x10-^* 


63 


18.18 


3 


1.37 


0.1 





68 


11.89 


10 


50 


-1.5614 


2x10-^* 


72 


44.19 


4 


3.05 


0.1 





86 


49.76 


10 


100 


-1.2812 


3x10-^* 


67 


58.38 


4 


3.51 


0.1 





178 


664.41 


20 


100 


13.5225 


2xKr* 


72 


56.68 


3 


2.83 


0.1 





217 


747.50 


20 


200 


-2.3946 


2xKr* 


76 


90.29 


4 


4.47 


0.1 





400 


17328.78 


30 


200 


3.3814 


3xKr* 


68 


81.06 


4 


6.62 


0.1 





290 


12937.88 


30 


300 


7.0023 


3xKr'* 


96 


203.63 


4 


10.47 


0.1 





# 


# 


40 


400 


17.2517 


3xKr'* 


77 


274.58 


4 


15.66 


0.1 





# 


# 


50 


500 


44.5623 


2xKr** 


100 


612.48 


4 


25.77 


0.1 





# 


# 



We first test Algorithms 1 and 2 with some random data in [—10, 10] and compare the 
performance with the built-in Matlab solver FMINCON for 10 problems. The data is created 
as follows: 

- Generate a random matrix M and compute Q := M^M + 0.5I„, where /„ is the identity 
matrix in R"^". 

- Vectors q, p, b and matrix A are random in [— 10, 10], and a = 10. 

- Generate a random matrix P^ in [— 10, 10] and then compute P := 0.5(Pr +^i^)- 

- The lower bound vector / and the upper bound vector u are given by (— 5, . . . , — 5)^ and 
(10, ... , 10)^, respectively. 

Since every symmetric matrix P can be decomposed as P = Pi — P2, where Pi and P2 are 
symmetric positive semidefinite (using spectral decompositions). The constraint 

-x^Px+p^x < a 
is expressed as a DC constraint: 

-x^Plx+p^x x^Pox < a, 

2 2 

where Pi = VE+V'^ and P2 = VZ-V^ with E+ = diag(a+) and £_ = diag(CT,."), CT;+ = 
max{CT,,0}, CT,.^ = max{— CT,',0}, CT/ is the /* eigenvalue of matrix P, and y is a matrix 
whose columns are formed by the eigenvectors of P. 

We implement Algorithms [T] and |2] in Matlab 7.8.0 (R2009a) running on a PC desk- 
top with Intel(R) Core(TM)2 Quad CPU Q6600 2.4GHz, 3Gb RAM. We use the same DC 
decomposition of the DC constraint in both algorithms. To solve the convex quadratic sub- 
problems, we use the CVX package (with Sedumi as a solverlJ. The tolerance is given by 
e = 10^^ and the penalty parameter fik is fixed to a certain value in Algorithm[2](see Tables 
[T]and[2]i. The computational results are reported in Table|T] For comparison, we solve three 
problems taken from I7j. The two first problems (Pi, P2) are in Chapter 3[7J[test problems 
1 and 2, respectively], while the last one (P3) is a VLSI design problem in Chapter 3|7][test 
problem 2]. The best-known solutions and optimal values of Pi, P2 are given in (TJ: 

4 = (579.31, 1359.97,5109.97, 182.02,295.6,217.98,286.42,395.60)^, /;= 7049.25, 
x*2 = (78,33,29.9953,45,36.7758)^, /I = 30665.5387, 



Available at: |http://cvxr.cQm/cvx/ 1 
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Table 2 Computational results of Algorithms ^and|2]for three nonconvex QP problems in 0. 



Problem Information 


AlgorithmtlJ(or|2j | 


N» 


[n.mi .ni2,ni3] 


otype 


r(in|7|) 


iter 


time 


W 


P 


/* 


Pi 


[8,0,3,3] 


In 


7049.25 


176 


74.36 


100 


1 X 10-^ 


7049.25352 


P2 


[5,0,0,6] 


nq 


-30665.5387 


5 


2.94 


810 


1 X 10-^' 


-30665.53892 


P3 


[12,5,2,6] 


nq 


146.25 


5 


2.76 








146.25000 



respectively. The optimal value of P3 is f^ = 146.25. Our computational results for these 
problems are reported in Table [21 which closely approximate to the best-known solution 
reported in fT\. 

The notations in Tables|T]and[2]include: n, m\, 1712, m^ are the size of the problems (vari- 
ables, linear equality, linear inequality and DC constraints, respectively), /* is the optimal 
value, otype is the type of the objective function (In is linear, nq is nonconvex quadratic), 
error is the quantity ||jr*^' "^'^ll' iter is the number of iterations, and time is the CPU 
time in seconds; jU and p are the penalty and the regularization parameters in Algorithm |2l 
respectively. The symbol # indicates that FMINCON exceeds the limit time Jniax = 4 hours. 



6.2 Example 2 



This example illustrates an application of Algorithmic] to solve mathematical programs with 
complementarity constraints presented in Section[2l 



' min f{x,y) 

s.t. Ax + By < a, 

x'^(Cx + Dy + e)=0, Cx + Dy + e>0, x>0, 
X £ Qf, y e Qy, 



(IVIPCC) 



where x 6 R"' is decision variable, y 6 R"' is parameter, / is convex with respect to x and 
y. A, B, C and D are given consistent matrices, a and e are given consistent vectors, and I2j, 
Qy are two convex sets in R"' and R"' , respectively. 

As in Example 12.21 of Section |2] we use a slack variable z for Cx + Dy + e > 0, the 
complementarity condition of dlVIFCCb is expressed equivalently to 



x>0,z>0,x'z<0, Cx + Dy-z + e = 0. 



(27) 



Let us define a new variable w 

|2 



{x^y, 



,T\T 



6 R"»' with «„. = Irix + riy, and denote by 
u[w) := ||(.ic,y,z)|p and v(m') := ||jc — z|p + ||y||^, the third condition of Ml\ is equivalent to 
a DC constraint u{w) — v{w) < 0. Note that u is strongly convex with parameter p" = 2 and 
V is only convex (not strongly convex). We also define 



^u := 



(jc^,/,z^)^eR"" \ Ax + By < a, Cx + Dy-z + e-- 
X e Q.X, y e 12,., Jc > 0, z > 



:0. 



(28) 



Since Q^ and Qy are convex, and the remaining constraints are linear, 12,^ is convex in R"" . 
Problem JlVIPCCt is reformulated as 



nun U{w):=f{x,y) 
s.t. u{w) — v(w) < 0, 



(29) 
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Table 3 Computational results of Algorithm |2]f or jMPCC) . 



N" 


Problem Information 


Algorithmic) | 




[m, n, 1] 


;t" 


/* 


error 


f easgap 


iter 


time 


P7 


(2,2,6) 


(40, 40) 


64.999 


7x 10-' 


5x10-" 


9 


9.91 


P? 


(2,2,2) 


(0, 0) 


7.095 X 10-'- 


1 X K)-^ 


4x10-'^ 


18 


13.00 


P? 


- 


(10, 0) 


1.351x10-" 


2 X K)-^ 


1 X 10- "' 


18 


12.82 


Pg 


- 


(5,5) 


1.294 X 10-" 


2x10-^ 


1 X 10-'" 


18 


12.83 


Pg 


- 


(0, 10) 


1.229 X 10-" 


2 X 10-^ 


Ix 10-" 


18 


12.91 


Pg 


- 


(10, 10) 


2.597 X 10-" 


8x W-" 


Ix K)-" 


19 


13.56 


Pio 


(4,4,12) 


(5,5,15,15) 


-6600 


3 X 10-=" 


3 X K)-** 


17 


15.76 



which coincides with ^. 

In this example, we implement Algorithmic] for solving three problems P7, P9 and Pio in 
l6lrproblems 7, 9 an d 10, resp ectively]. The parameter ^k is fixed to jU = IQ-' . To solve the 
convex subproblems P(j«^;/i) we also use the CVX package with the Sedumi solver. For a 
given tolerance e = IQ-'', the computational results are presented in Table [3] which closely 
approximate to the results given in f^. The solutions reported by Algorithm |2] for P7, P9 and 
Pio are 

4^ = (25.00125,30.00000)^, 4, = (10,5)^ 
and4j^ = (7.515728,3.77360,11.48427,17.22640)^, 

respectively. Algorithm [T] failed in this case because the set of interior points intD of the 
feasible set D is empty. 



7 Conclusion 



The main aim of this paper is to investigate the relation between sequential convex program- 
ming (SCP) I18II24I and DC programming 1 1 , 2 23 1. We have provided a variant of the SCP 
algorithm for finding local minimizers of a nonconvex programming problem with DC con- 
straints. We have proved a global convergence theorem for this particular algorithm. Then 
we have addressed some extensions and proposed a relaxation technique to handle possibly 
inconsistent linearizations. Although finding a DC decomposition of a certain DC function is 
in general still a hard problem, in some applications (as we have shown in the examples) it is 
available or easy to compute. We have not concentrated on the local convergence. However, 
under mild assumptions, it had been proved in [24] that the SCP method converges linearly 
to a KKT point of the original problem. Applications to nonconvex quadratic programming 
problems as well as mathematical programming problems with complementarity constraints 
have been presented through two numerical examples. 
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